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This paper presents algorithms for hierarchical, agglomerative clustering which 
CN perform most efficiently in the general-purpose setup that is given in modern 

O | standard software. Requirements are: (1) the input data is given by pairwise dis- 

D similarities between data points, but extensions to vector data are also discussed 

(2) the output is a "stepwise dendrogram", a data structure which is shared by 
£S) all implementations in current standard software. We present algorithms (old and 

t-H new) which perform clustering in this setting efficiently, both in an asymptotic 

worst-case analysis and from a practical point of view. The main contributions of 
this paper are: (1) We present a new algorithm which is suitable for any distance 
update scheme and performs significantly better than the existing algorithms. (2) 
We prove the correctness of two algorithms by Rohlf and Murtagh, which is neces- 
sary in each case for different reasons. (3) We give well-founded recommendations 
for the best current algorithms for the various agglomerative clustering schemes. 
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1^ 1 Introduction 

m 

Hierarchical, agglomerative clustering is an important and well-established technique in un- 
Q^ supervised machine learning. Agglomerative clustering schemes start from the partition of 

the data set into singleton nodes and merge step by step the current pair of mutually closest 
nodes into a new node until there is one final node left, which comprises the entire data set. 
Various clustering schemes share this procedure as a common definition, but differ in the way 
in which the measure of inter-cluster dissimilarity is updated after each step. The seven most 
common methods are termed single, complete, average (UPGMA), weighted (WPGMA, Mc- 
Quitty), Ward, centroid (UPGMC) and median (WPGMC) linkage (see Everitt et al., 2011, 
Table 4.1). They are implemented in standard numerical and statistical software such as R 
(R Development Core Team, 2011), MATLAB (The MathWorks, Inc., 2011), Mathematica 
(Wolfram Research, Inc., 2010), SciPy (Jones ct al., 2001). 

The stepwise, procedural definition of these clustering methods directly gives a valid but 
inefficient clustering algorithm. Starting with Gower's and Ross's observation (Gower and 
Ross, 1969) that single linkage clustering is related to the minimum spanning tree of a graph 
in 1969, several authors have contributed algorithms to reduce the computational complexity 
of agglomerative clustering, in particular Sibson (1973), Rohlf (1973), Anderberg (1973, page 
135), Murtagh (1984), Day and Edelsbrunner (1984, Table 5). 
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Even when software packages do not use the inefficient primitive algorithm (as SciPy (Eads, 
2007) and the R (R Development Core Team, 2011) methods hclust and agnes do), the 
author found that implementations largely use suboptimal algorithms rather than improved 
methods suggested in theoretical work. This paper is to advance the theory further up to a 
point where the algorithms can be readily used in the standard setting, and in this way bridge 
the gap between the theoretical advances that have been made and the existing software 
implementations, which are widely used in science and industry. 

The main contributions of this paper are: 

• We present a new algorithm which is suitable for any distance update scheme and per- 
forms significantly better than the existing algorithms for the "centroid" and "median" 
clustering schemes. 

• We prove the correctness of two algorithms, a single linkage algorithm by Rohlf (1973) 
and Murtagh's nearest-neighbor-chain algorithm (Murtagh, 1985, page 86). These 
proofs were still missing, and we detail why the two proofs are necessary, each for 
different reasons. 

• These three algorithms (together with an alternative by Sibson, 1973) are the best 
currently available ones, each for its own subset of agglomerative clustering schemes. 
We justify this carefully, discussing potential alternatives. 

The specific class of clustering algorithms which is dealt with in this paper has been char- 
acterized by the acronym SAHN (sequential, agglomerative, hierarchic, nonoverlapping meth- 
ods) by Sneath and Sokal (1973, §5.4, 5.5). The procedural definition (which is given in 
Figure 1 below) is not the only possibility for a SAHN method, but this method together with 
the seven common distance update schemes listed above is most widely used. The scope of 
this paper is contained further by practical considerations: We consider methods here which 
comply to the input and output requirements of the general-purpose clustering functions in 
modern standard software: 

• The input to the algorithm is the list of (^) pairwise dissimilarities between N points. 
(We mention extensions to vector data in Section 6.) 

• The output is a so called stepwise dendrogram (see Section 2.2), in contrast to laxly 
specified output structure or weaker notions of (non-stepwise) dendrograms in earlier 
literature. 

The first item has always been a distinctive characteristic to previous authors since the input 
format broadly divides into the stored matrix approach (Anderberg, 1973, §6.2) and the stored 
data approach (Anderberg, 1973, §6.3). In contrast, the second condition has not been given 
attention yet, but we will see that it affects the validity of algorithms. 

We do not aim to present and compare all available clustering algorithms but build upon 
the existing knowledge and present only the algorithms which we found the best for the given 
purpose. For reviews and overviews we refer to Rohlf (1982), Murtagh (1983, 1985), Gordon 
(1987, §3.1), Jain and Dubes (1988, §3.2), Day (1996, §4.2), Hansen and Jaumard (1997). 
Those facts about alternative algorithms which are necessary to complete the discussion and 
which are not covered in existing reviews are collected in Section 5. 

The paper is structured as follows: 

Section 2 contains the definitions for input and output data structures as well as specifica- 
tions of the distance update formulas and the "primitive" clustering algorithm. 
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Section 3 is the main section of this paper. We present and discuss three algorithms: our 
own "generic algorithm", Murtagh's nearest-neighbor-chain algorithm and Rohlf's algorithm 
based on the minimum spanning tree of a graph. We prove the correctness of these algorithms. 

Section 4 discusses the complexity of the algorithms, both as theoretical, asymptotic com- 
plexity in Section 4.1 and by use-case performance experiments in Section 4.2. We conclude 
this section by recommendations on which algorithm is the best one for each distance update 
scheme, based on the preceding analysis. 

Section 5 discusses alternative algorithms, and Section 6 gives a short outlook on extending 
the context of this paper to vector data instead of dissimilarity input. The paper ends with 
a brief conclusion in Section 7. 

The algorithms in this paper have been implemented in C++ by the author and are available 
with interfaces to the statistical software R and the programming language Python (van 
Rossum et al.). This implementation is presented elsewhere (Milliner, 2011). 

2 Data structures and the algorithmic definition of 
SAHN clustering methods 

In this section, we recall the common algorithmic (procedural) definition of the SAHN clus- 
tering methods which demarcate the scope of this paper. Before we do so, we concretize the 
setting further by specifying the input and output data structures for the clustering methods. 
Especially the output data structure has not been specifically considered in earlier works, but 
nowadays there is a de facto standard given by the shared conventions in the most widely 
used software. Hence, we adopt the setting from practice and specialize our theoretical con- 
sideration to the modern standard of the stepwise dendrogram. Later, Section 5 contains an 
example of how the choice of the output data structure affects the result which algorithms 
are suitable and/or most efficient. 

2.1 Input data structure 

The input to the hierarchical clustering algorithms in this paper is always a finite set together 
with a dissimilarity index (sec Hansen and Jaumard, 1997, § 2.1). 

Definition. A dissimilarity index on a set S is a map d : S x S — > [0, oo) which is reflexive 
and symmetric, i.e. we have d(x, x) = and d(x, y) = d(y, x) for all x,y G S. 

A metric on S is certainly a dissimilarity index. In the scope of this paper, we call the 
values of d distances in a synonymous manner to dissimilarities, even though they are not 
required to fulfill the triangle inequalities, and dissimilarities between different elements may 
be zero. 

If the set S has N elements, a dissimilarity index is given by the (^) pairwise dissimilarities. 
Hence, the input size to the clustering algorithms is Q(N 2 ). Once the primitive clustering 
algorithm is specified in Section 2.4, it is easy to see that the hierarchical clustering schemes 
are sensitive to each input value. More precisely, for every input size N and for every index 
pair i ^ j, there are two dissimilarities which differ only at position (i,j) and which produce 
different output. Hence, all input values must be processed by a clustering algorithm, and 
therefore the run-time is bounded below by £l(N 2 ). 

This bound applies to the general setting when the input is a dissimilarity index. In a 
different setting, the input could also be given as N points in a normed vector space of 
dimension D (the "stored data approach", Anderberg, 1973, §6.3). This results in an input 
size of Q(ND), so that the lower bound does not apply for clustering of vector data. See 
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Section 6 for a discussion to which extent the algorithms in this paper can be used in the 
"stored data approach". 

2.2 Output data structures 

The output of a hierarchical clustering procedure is traditionally a dendrogram. The term 
"dendrogram" has been used with three different meanings: a mathematical object, a data 
structure and a graphical representation of the former two. In the course of this section, we 
define a data structure and call it stepwise dendrogram. A graphical representation may be 
drawn from the data in one of several existing fashions. The graphical representation might 
lose information (e.g. when two merges happen at the same dissimilarity value), and at the 
same time contain extra information which is not contained in the data itself (like a linear 
order of the leaves). 

In the older literature, e.g. Sibson (1973), a dendrogram (this time, as a mathematical 
object) is rigorously defined as a piecewise constant, right-continuous map D: [0, oo) — > V(S), 
where V(S) denotes the partitions of S, such that 

• D(s) is always coarser than or equal to D(t) for s > t, 

• D(s) eventually becomes the one-set partition {S} for large s. 

A dendrogram in this sense with the additional condition that D(0) is the singleton partition 
is in one-to-one correspondence to an ultrametric on S (Johnson, 1967, § I). The ultrametric 
distance between x and y is given by n(x,y) — min{s > | x ~ y in D(s)}. Conversely, 
the partition at level s > in the dendrogram is given by the equivalence relation 

y) < s - Sibson's "pointer representation" and "packed representation" (Sibson, 1973, § 4) 
are examples of data structures which allow the compact representation of a dendrogram or 
ultrametric. 

In the current software, however, the output of a hierarchical clustering procedure is a 
different data structure which conveys potentially more information. We call this a stepwise 
dendrogram. 

Definition. Given a finite set So with cardinality N = \Sq\, a stepwise dendrogram is a list 
of N — 1 triples (o», 6», 5j) (i = 0, . . . , N — 2) such that Si € [0, oo) and a,, bi G Si, where Si+i 
is recursively defined as (Si \ {a,, bi}) U n* and n, ^ S\ {a i7 bi} is a label for a new node. 

This has the following interpretation: The set So are the initial data points. In each step, 
n, is the new node which is formed by joining the nodes a, and bi at the distance Si. The 
order of the nodes within each pair (aj,6j) does not matter. The procedure contains N — 1 
steps, so that the final state is a single node which contains all N initial nodes. 

(The mathematical object behind this data structure is a sequence of N+ 1 distinct, nested 
partitions from the singleton partition to the one-set partition, together with a nonnegative 
real number for each partition. We do not need this abstract point of view here, though.) 

The identity of the new labels rii is not part of the data structure; instead it is assumed that 
they are generated according to some rule which is part of the specific data format convention. 
In view of this, it is customary to label the initial data points and the new nodes by integers. 
For example, the following schemes are used in software: 

• R convention: So ■= (— 1, • • • , — N), new nodes: (1, . . . ,N — 1) 

• SciPy convention: S (0, . . . , N - 1), new nodes: (N, . . . , 2N — 2) 

• MATLAB convention: S Q :=(!,..., N), new nodes: (N + 1, . . . , 2N - 1) 
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We regard a stepwise dendrogram both as an ((N — 1) x 3)-matrix or as a list of triples, 
whichever is more convenient in a given situation. 

If the sequence (Si) in Section 2.2 is non-decreasing, one says that the stepwise dendrogram 
does not contain inversions, otherwise it does. 

In contrast to the first notion of a dendrogram above, a stepwise dendrogram can take 
inversions into account, which an ordinary dendrogram cannot. Moreover, if more than two 
nodes are joined at the same distance, the order of the merging steps does matter in a stepwise 
dendrogram. Consider e.g. the following data sets with three points: 



(A) x 

2.0 / \ 2.0 
/ \ 
•i ^» 

X\ 3.0 x 2 



(B) 

2.0 / \ 2.0 
/ \ 

^» 

X 2 3.0 x 



(C) 

2.0 / \ 2.0 
/ \ 
^» 

X 3.0 X\ 



The array 

"0, 1, 2.0" 

2, 3, 2.0 

is a valid output (SciPy conventions) for single linkage clustering on the data sets (A) and 
(B) but not for (C). Even more, there is no stepwise dendrogram which is valid for all three 
data sets simultaneously. On the other hand, the non-stepwise single linkage dendrogram is 
the same in all cases: 

m fW.M.M if*<2 pictorially: p-^" 

Y{xq,x\,X2} if s > 2 x xi x 2 

Hence, a stepwise dendrogram conveys slightly more information than a non-stepwise dendro- 
gram in the case of ties (i.e. when more than one merging step occurs at a certain distance). 
This must be taken into account when we check the correctness of the algorithms. Although 
this complicates the proofs in Sections 3.3 and 3.2 and takes away from the simplicity of the 
underlying ideas, it is not a matter of hairsplitting: E.g. Sibson's SLINK algorithm (Sibson, 
1973) for single linkage clustering works flawlessly if all distances are distinct but produces 
the same output on all data sets (A), (B) and (C). Hence, the output cannot be converted 
into a stepwise dendrogram. See Section 5 for further details. 



2.3 Node labels 

The node labels rij in a stepwise dendrogram may be chosen as unique integers according to 
one of the schemes described in the last section. In an implementation, when the dissimilarities 
are stored in a large array in memory, it is preferable if each node label rii for the joined cluster 
reuses one of the indices ai,bi of its constituents, so that the dissimilarities can be updated 
in-place. Since the clusters after each row in the dendrogram form a partition of the initial set 
Sq, we can identify each cluster not only by its label but also by one of its members. Hence, if 
the new node label rij is chosen among a.i,bi, this is sufficient to reconstruct the partition at 
every stage of the clustering process, and labels can be converted to any other convention in 
a postprocessing step. Generating unique labels from cluster representatives takes only Q(N) 
time and memory with a suitable union-find data structure. See Section 3.2 and Figure 5 for 
details. 



5 



Figure 1 Algorithmic definition of a hierarchical clustering scheme. 



l: procedure Primitive_clustering(S', d) 

2: N*-\S\ 

3: L<-[] 

4: size[x] 1 for all x <G S 

5: for i <— 0, . . . , N — 2 do 
6: (a, 6) «- argmin (SxS)XA d 

7: Append (a, &, d[a, b]) to i. 

8: S 4- S \ {a, 6} 

9: Create a new node label n S. 

10: Update d with the information 



> S: node labels, d: pairwise dissimilarities 
> Number of input nodes 
> Output list 



c/[n,x] = d[a;,n] = FORMULA(d[a, x], d[b, x], d[a, b], size[a], size[b], size[x]) 

for all x £ S. 
11: size[n] size[a] + size[b] 

12: S^SU{n} 
13: end for 

14: return L > the stepwise dendrogram, an ((N — 1) x 3)-matrix 

15: end procedure 

(As usual, A denotes the diagonal in the Cartesian product 5x5.) 



2.4 The primitive clustering algorithm 

The solution that we expect from a hierarchical clustering algorithm is defined procedurally. 
All algorithms in this paper are measured against the primitive algorithm in Figure 1. We state 
it in a detailed form to point out exactly which information about the clusters is maintained: 
the pairwise dissimilarities and the number of elements in each cluster. 

The function Formula in line 10 is the distance update formula, which returns the distance 
from a node x to the newly formed node a U b in terms of the dissimilarities between clusters 
a, b and x and their sizes. The table in Figure 2 lists the formulas for the common distance 
update methods. 

For five of the seven formulas, the distance between clusters does not depend on the order 
which the clusters were formed by merging. In this case, we also state closed, non-iterative 
formulas for the cluster dissimilarities in the third row in Figure 2. The distances in the 
"weighted" and the "median" update scheme depend on the order, so we cannot give non- 
iterative formulas. 

The "centroid" and "median" formulas can produce inversions in the stepwise dendrograms; 
the other five methods cannot. This can be checked easily: The sequence of dissimilarities at 
which clusters are merged in Figure 1 cannot decrease if the following condition is fulfilled for 
all disjoint subsets I,J,KcSo- 

d(I,J)<mm{d(I,K),d(J,K)} => d(I, J) < d(I U J, K) 

On the other hand, configurations with inversion in the "centroid" and "median" schemes 
can be easily produced, e.g. three points near the vertices of an equilateral triangle in M 2 . 

The primitive algorithm takes ©(TV 3 ) time since in the i-th iteration of N — 1 in total, all 
( N ~2~ l ) *= ©(* 2 ) pairwise distances between the N — i nodes in S are searched. 

Note that the stepwise dendrogram from a clustering problem (5, d) is not always uniquely 
defined, since the minimum in line 6 of the algorithm might be attained for several index 
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Figure 2 Agglomerative clustering schemes. 



Name 



Distance update formula 
Formula for d(I U J, K) 



Cluster dissimilarity 
between clusters A and B 



single 
complete 

average 

weighted 

Ward 

centroid 
median 



mm(d(I,K),d(J,K)) 

max(d(I,K),d(J,K)) 

md(I, K) + njd(J, K) 
m + nj 

d(I,K)+d(J,K) 



(nj + n K )d(I, K) + (nj + n K )d(J, K) - n K d(I, J) 
ni + nj + n K 

njd(I, K) + njd(J, K) njnjd{I,J) 



m + nj 



(rn + nj) 2 



d(I,K) d{J,K) _ d(I,J) 



min d\a,b] 

aGA,b<£B 

max d\a,b] 

a£A,b£B 



\A\\B 



a£A beB 



2\A\\B\ 
\A\ + \B\ 



\\CA - C B 2 



\CA - C B \\2 



\WA - W B \\2 



Legend: Let /, J be two clusters joined into a new cluster, and let K be any other cluster. 
Denote by nj, nj and uk the sizes of (i.e. number of elements in) clusters /, J, K, respectively. 

The update formulas for the "Ward", "centroid" and "median" methods assume that the input 
points are given as vectors in Euclidean space with the Euclidean distance as dissimilarity 
measure. The expression cx denotes the centroid of a cluster X. The point wx is defined 
iteratively and depends on the clustering steps: If the cluster L is formed by joining I and J, 
we define Wl as the midpoint \(wi + wj). 

All these formulas can be subsumed (for squared Euclidean distances in the three latter cases) 
under a single formula 

d(I U J, K) := ai d(I, K) + ajd(J, K) + (3d(I, J) + j\d(I, K) - d(J, K)\, 

where the coefficients ai,aj,j3 may depend on the number of elements in the clusters /, J 
and K. For example, aj = otj = |,/3 = 0, 7 = — \ gives the single linkage formula. All 
clustering methods which use this formula are combined under the name "flexible" in this 
paper, as introduced by Lance and Williams (1967). 

References: Lance and Williams (1967), Kaufman and Rousseeuw (1990, §5.5.1) 
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pairs. We consider every possible output of Primitive_CLUSTERING under any choices of 
minima as a valid output. 

3 Algorithms 

In the main section of this paper, we present three algorithms which are the most efficient ones 
for the task of SAHN clustering with the stored matrix approach. Two of the algorithms were 
described previously: The nearest-neighbor chain ("NN-chain") algorithm by Murtagh (1985, 
page 86), and an algorithm by Rohlf (1973), which we call "MST-algorithm" here since it is 
based on Prim's algorithm for the minimum spanning tree of a graph. Both algorithms were 
presented by the respective authors, but for different reasons each one still lacks a correctness 
proof. Sections 3.2 and 3.3 state the algorithms in a way which is suitable for modern standard 
input and output structures and supply the proofs of correctness. 

The third algorithm in Section 3.1 is a new development based on Anderberg's idea to 
maintain a list of nearest neighbors for each node Anderberg (1973, pages 135-136). While 
we do not show that the worst-case behavior of our algorithm is better than the 0(N 3 ) worst- 
case complexity of Anderberg's algorithm, the new algorithm is for all inputs at least equally 
fast, and we show by experiments in Section 4.2 that the new algorithm is considerably faster 
in practice since it cures Anderberg's algorithm from its worst-case behavior at random input. 

As we saw in the last section, the solution to a hierarchical clustering task does not have a 
simple, self-contained specification but is defined as the outcome of the "primitive" clustering 
algorithm. The situation is complicated by the fact that the primitive clustering algorithm 
itself is not completely specified: if a minimum inside the algorithm is attained at more than 
one place, a choice must be made. We do not require that ties are broken in a specific 
way; instead we consider any output of the primitive algorithm under any choices as a valid 
solution. Each of the "advanced" algorithms is considered correct if it always returns one of 
the possible outputs of the primitive algorithm. 

3.1 The generic clustering algorithm 

The most generally usable algorithm is described in this section. We call it Generic_ 
linkage since it can be used with any distance update formula. It is the only algorithm among 
the three in this paper which can deal with inversions in the dendrogram. Consequentially, 
the "centroid" and "median" methods must use this algorithm. 

The algorithm is presented in Figure 3. It is a sophistication of the primitive clustering 
algorithm and of Anderberg's approach (Anderberg, 1973, pages 135-136). Briefly, candidates 
for nearest neighbors of clusters are cached in a priority queue to speed up the repeated 
minimum searches in line 6 of Primitive_CLUSTERING. 

For the pseudocode in Figure 3, we assume that the set S are integer indices from to 
N — 1. This is the way in which it may be done in an implementation, and it makes the 
description easier than for an abstract index set S. In particular, we rely on an order of the 
index set (see e.g. line 6: the index ranges over all y > x). 

There are two levels of sophistication from the primitive clustering algorithm to our generic 
clustering algorithm. In a first step, one can maintain a list of nearest neighbors for each 
cluster. For the sake of speed, it is enough to search for the nearest neighbors of a node x 
only among the nodes with higher index y > x. Since the dissimilarity index is symmetric, this 
list will still contain a pair of closest nodes. The list of nearest neighbors speeds up the global 
minimum search in the i-th step from ( JV ~ 2 * _1 ) comparisons to N — i — 1 comparisons at the 
beginning of each iteration. However, the list of nearest neighbors must be maintained: if the 
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Figure 3 The generic clustering algorithm. 

1: procedure Generic linkage(]V, d) t> N: input size, d: pairwise dissimilarities 

2: S<- (0,...,JV- 1) 

3: L 4- [] > Output list 

4: size [a;] 4- 1 for all x £ S 

5: for x in S\ {N — 1} do > Generate the list of nearest neighbors. 

6: n_nghbr[x] 4- argmin J/>;r d[x, y] 

7: mindist[x] 4- d[x, n_nghbr[x\] 

8: end for 

9: Q 4- (priority queue of indices in S \ {N — 1}, keys are in mindist) 

10: for i 4- 1, . . . , N - 1 do > Main loop. 

11: a 4— (minimal element of Q) 

12: & n_ng/i&r[a] 

13: 5 <— mmdisi[a] 

14: while 5 ^ d{a, b] do > Recalculation of nearest neighbors, if necessary. 

15: n_nghbr[a] argmin 2;>a d[a, x] 

16: Update mindist and Q with (a, d[a, n_nghbr[a]]) 

17: a (minimal element of Q) 

18: <— n_nghbr[a] 

19: <5 <— ™W£s£[a] 

20: end while 

21: Remove the minimal clement a from Q. 

22: Append (a, b, S) to L. > Merge the pairs of nearest nodes. 

23: size[6] 4- size[a] + size[b) > Re-use b as the index for the new node. 

24: 5'^S'\{a} 

25: for x in S \ {b} do > Update the distance matrix. 

26: d[x, b] d[b, x] FORMULA(d[a, x], d[b, x],d[a, b], size[a], size [b] , size [x] ) 

27: end for 

28: for x in 5 such that x < a do > Update candidates for nearest neighbors. 

29: if n_nghbr[x] = a then > Deferred search; no nearest 

30: n_nghbr[x] b > neighbors are searched here. 

31: end if 

32: end for 

33: for x in S such that x < b do 

34: if d[x, b] < mindist[x] then 

35: n_nghbr[x] 6 

36: Update mindist and Q with (x,d[x,b\) t> Preserve a lower bound. 

37: end if 

38: end for 

39: n_nghbr[b] <— argmin a , >b d[b, x] 

40: Update mindist and Q with (6, d[b, n_nghbr[b]]) 

41: end for 

42: return L > The stepwise dendrogram, an ((TV — 1) x 3)-matrix. 
43: end procedure 
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nearest neighbor of a node x is one of the clusters a, b which are joined, then it is sometimes 
necessary to search again for the nearest neighbor of x among all nodes y > x. Altogether, 
this reduces the best-case complexity of the clustering algorithm from ©(TV 3 ) to 0(iV 2 ), while 
the worst case complexity remains 0(N 3 ). This is the method that Anderberg suggested in 
(Andcrberg, 1973, pages 135-136). 

On a second level, one can try to avoid or delay the nearest neighbor searches as long as 
possible. Here is what the algorithm Generic_linkage does: It maintains a list n_nghbr 
of candidates for nearest neighbors, together with a list mindist of lower bounds for the 

distance to the true nearest neighbor. If the distance d[x, n nghbr[x]] is equal to mindist[x], 

we know that we have the true nearest neighbor, since we found a realization of the lower 
bound; otherwise the algorithm must search for the nearest neighbor of x again. 

To further speed up the minimum searches, we also make the array mindist into a priority 
queue, so that the current minimum can be found quickly. We require a priority queue Q with 
a minimal set of methods as in the list below. This can be implemented conveniently by a 
binary heap (see Cormen et al., 2009, Chapter 6). We state the complexity of each operation 
by the complexity for a binary heap. 

• Queue(w): Generate a new queue from a vector v of length \v\ — N. Return: an object 
Q. Complexity: O(N). 

• Q.Argmin: Return the index to a minimal value of v. Complexity: O(l). 

• Q. Remove Min: Remove the minimal element from the queue. Complexity: O(logiV). 

• Q.UPDATE(i, x): Assign v[i] x and update the queue accordingly. Complexity: 
O (log N). 

We can now describe the Generic LINKAGE algorithm step by step: Lines 5 to 8 search 

the nearest neighbor and the closest distance for each point x among all points y > x. This 
takes 0(N 2 ) time. In line 9, we generate a priority queue from the list of nearest neighbors 
and minimal distances. 

The main loop is from line 10 to the end of the algorithm. In each step, the list L for 
a stepwise dendrogram is extended by one row, in the same way as the primitive clustering 
algorithm does. 

Lines 11 to 20 find a current pair of closest nodes. A candidate for this is the minimal index 
in the queue (assigned to a), and its candidate for the nearest neighbor b := n_nghbr[a\. If 
the lower bound mindist[a] is equal to the actual dissimilarity d[a, b], then we are sure that 
we have a pair of closest nodes and their distance. Otherwise, the candidates for the nearest 
neighbor and the minimal distance are not the true values, and we find the true values for 

n nghbr[a] and mindist[a] in line 15 in O(N) time. We repeat this process and extract the 

minimum among all lower bounds from the queue until we find a valid minimal entry in the 
queue and therefore the actual closest pair of points. 

This procedure is the performance bottleneck of the algorithm. The algorithm might be 
forced to update the nearest neighbor O(N) times with an effort of O(N) for each of the 
O(N) iterations, so the worst-case performance is bounded by 0(N 3 ). In practice, the inner 
loop from lines 14 to 20 is executed less often, which results in faster performance. 

Lines 22 to 27 are nearly the same as in the primitive clustering algorithm. The only 
difference is that we specialize from an arbitrary label for the new node to re-using the index 
b for the joined node. The index a becomes invalid, and we replace any nearest-neighbor 
reference to a by a reference to the new cluster b in lines 28 to 32. Note that the array 
n_nghbr contains only candidates for the nearest neighbors, so we could have written any 



10 



valid index here; however, for the single linkage method, it makes sense to choose b: if the 
nearest neighbor of a node was at index a, it is now at b, which represents the join a U b. 

The remaining code in the main loop ensures that the array n nghbr still contains lower 

bounds on the distances to the nearest neighbors. If the distance from the new cluster x to 
a cluster b < x is smaller than the old bound for b, we record the new smallest distance and 
the new nearest neighbor in lines 34 to 37. 

Lines 39 and 40 finally find the nearest neighbor of the new cluster and record it in the 
arrays n_nghbr and mindist and the queue Q 

The main idea behind this approach is that invalidated nearest neighbors are not re- 
computed immediately. Suppose that the nearest neighbor of a node x is far away from 
x compared to the global closed pair of nodes. Then it does not matter that we do not know 
the nearest neighbor of x, as long as we have a lower bound on the distance to the nearest 
neighbor. The candidate for the nearest neighbor might remain invalid, and the true distance 
might remain unknown for many iterations, until the lower bound for the nearest-neighbor 
distance has reached the top of the queue Q. By then, the set of nodes S might be much 
smaller since many of them were already merged, and the algorithm might have avoided many 
unnecessary repeated nearest-neighbor searches for x in the meantime. 

This concludes the discussion of our generic clustering algorithm; for the performance see 
Section 4. Our explanation of how the minimum search is improved also proves the correctness 
of the algorithm: Indeed, in the same way as the primitive algorithm does, the Generic_ 
linkage algorithm finds a pair of globally closest nodes in each iteration. Hence the output 
is always the same as from the primitive algorithm (or more precisely: one of several valid 
possibilities if the closest pair of nodes is not unique in some iteration) . 

3.2 The nearest-neighbor-chain algorithm 

In this section, we present and prove correctness of the nearest-neighbor-chain algorithm 
(shortly: NN-chain algorithm), which was described by Murtagh (1985, page 86). This algo- 
rithm can be used for the "single", "complete", "average", "weighted" and "Ward" methods. 

The NN-chain algorithm is presented in Figure 4 as NN-CHAIN-linkage. It consists of the 
core algorithm NN-CHAIN-CORE and two postprocessing steps. Because of the postprocessing, 
we call the output of NN-CHAIN-CORE an unsorted dendrogram. The unsorted dendrogram 
must first be sorted row-wise, with the dissimilarities in the third column as the sorting key. 
In order to correctly deal with merging steps which happen at the same dissimilarity, it is 
crucial that a stable sorting algorithm is employed, i.e. one which preserves the relative order 
of elements with equal sorting keys. At this point, the first two columns of the output array 
L contain the label of a member of the respective cluster, but not the unique label of the 
node itself. The second postprocessing step is to generate correct node labels from cluster 
representatives. This can be done in 0(iV) time with a union-find data structure. Since this 
is a standard technique, we do not discuss it here but state an algorithm in Figure 5 for the 
sake of completeness. It generates integer node labels according to the convention in SciPy 
but can easily be adapted to follow any convention. 

We prove the correctness of the NN-chain algorithm in this paper for two reasons: 

• We make sure that the algorithm resolves ties correctly, which was not in the scope of 
earlier literature. 

• Murtagh claims (Murtagh, 1984, page 111), (Murtagh, 1985, bottom of page 86) that 
the NN-chain algorithm works for any distance update scheme which fulfills a certain 
"reducibility property" 
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d{I, J) < mm{d(I, K),d(J, K)} => mm{d(I,K),d(J,K)} < d(I U J,K) (1) 



for all disjoint nodes I,J,K at any stage of the clustering (Murtagh, 1984, §3), (Mur- 
tagh, 1985, § 3.5) . This is false. 1 We give a correct proof which also shows the limitations 
of the algorithm. In Murtagh's papers (Murtagh, 1984, 1985), it is not taken into account 
that the dissimilarity between clusters may depend on the order of clustering steps; on 
the other hand, it is explicitly said that the algorithm works for the "weighted" scheme, 
in which dissimilarities depend on the order of the steps. 

Since there is no published proof for the NN-chain algorithm but claims which go beyond 
what the algorithm can truly do, it is necessary to establish the correctness by a strict proof: 

Theorem 1. Fix a distance update formula. For any sequence of merging steps and any four 
disjoint clusters I, J, K, L resulting from these steps, require two properties from the distance 
update formula: 

• It fulfills the reducibility property (1). 

• The distance d(I U J, K U L) is independent of whether (I, J) are merged first and then 
(K, L) or the other way round. 

Then the algorithm NN-CHAIN-LINKAGE produces valid stepwise dendrograms for the given 
method. 

Proposition 2. The "single", "complete", "auerage", "weighted" and "Ward" distance update 
formulas fulfill the requirements of Theorem 1. 

Proof of Theorem 1. We prove the theorem by induction in the size of the input set S. The 
induction start is trivial since a dendrogram for a one-point set is empty. 

We call two nodes a,b e S reciprocal nearest neighbors ("pairwise nearest neighbors" in the 
terminology of Murtagh, 1985) if the distance d[a, b] is minimal among all distances from a to 
points in S, and also minimal among all distances from b: 



d[a, b] — min d[a, x] = min d[b, x] 

x^a x^b 



Every finite set S with at least two elements has at least one pair of reciprocal nearest 
neighbors, namely a pair which realizes the global minimum distance. 

The list chain is in the algorithm constructed in a way such that every element is a nearest 
neighbor of its predecessor. If chain ends in [. . . , b, a, b], we know that a and b are reciprocal 
nearest neighbors. The main idea behind the algorithm is that reciprocal nearest neighbors 
a,b always contribute a row (a, b, d[a, b]) to the stepwise dendrogram, even if they are not 
discovered in ascending order of dissimilarities. 



1 For example, consider the distance update formula d(I U J, K) := d(I, K) + d(J, K) + d(I, J). This formula 
fulfills the reducibility condition. Consider the following distance matrix between five points in the first 

column below. The Primitive clustering algorithm produces the correct stepwise dendrogram in the 

middle column. However, if the point A is chosen first in line 7 of NN-CHAIN-CORE, the algorithm outputs 
the incorrect dendrogram in the right column. 

(C, D, 1) (C, D, 1) 

(A, B, 3) (A, B, 3) 

(AB, CD, 27) (CD, E, 28) 

(ABCD, E, 85) {AB, CDE, 87) 
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Figure 4 The nearest-neighbor clustering algorithm. 



procedure NN-CHAIN-linkage(S', d) > S: node labels, d: pairwise dissimilarities 

L <- NN-CHAIN-CORE(iV, d) 

Stably sort L with respect to the third column. 

L <— Label(L) > Find node labels from cluster representatives, 

return L 
end procedure 

procedure NN-CHAIN-CORe(5', d) > S: node labels, d: pairwise dissimilarities 

S^(0,...,N-1) 
chain = [ ] 

size[x] 1 for all x £ S 
while \S\ > 1 do 

if length (c/iain) < 3 then 

a <r- (any element of S) > E.g. S[0] 

chain <— [a] 

b «- (any element of S \ {a}) > E.g. S[l] 

else 

a 4— chain[— 4] 
<— chain[— 3] 

Remove chain[— 1], chain[—2] and chain[— 3] > Cut the tail (x,y,x). 

end if 
repeat 

c <— argmin a .^ a a] with preference for 6 
a, 6 c, a 
Append a to ckm 

until length(c/iam) > 3 and a — chain[— 3] > a, 6 are reciprocal 

Append (a, 6, d[a, &]) to L > nearest neighbors. 

Remove a, 6 from 5 
ra <— (new node label) 
size[n] 4— size [a] + size [6] 
Update d with the information 

d[n, x] = rf[a;,n] = FORMULA(<i[a, x], d[b, x], d[a, b], size[a], size[b], size[x]) 

for all x £ S. 
S 4-SU{n} 
end while 

return L > an unsorted dendrogram 

end procedure 

(We use the Python index notation: chain[— 2] is the second-to-last element in the list chain.) 
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Figure 5 A union- find data structure suited for the output conversion. 



1: procedure Label(L) 

2: L'<-[] 

3: N <— (number of rows in L) + 1 > Number of initial nodes. 

4: U <- new Union-Find (iV) 

5: for (a, 6, S) in L do 

6: Append (£/.EFFiCiENT-FiND(a), L/.Efficient-Find(6), S) to L' 

7: {/.UNION(<Z, 6) 

8: end for 

9: return L' 
10: end procedure 

li: class Union-Find 

12: method CONSTRUCTOR(iV) > N is the number of data points. 

13: parent <— new int[27V — 1] 

14: parent[0, ...,2N-2] <- None 

15: nextlabel ^ N > SciPy convention: new labels start at AT 

16: end method 

17: method UNiON(m, n) 
18: parent [to] = nextlabel 

19: parent[n] = nextlabel 

20: nextlabel nextlabel +1 > SciPy convention: number new labels consecutively 

21: end method 

22: method FlND(n) > This works but the search process is not efficient. 

23: while parent [n] is not None do 

24: n <— parent[n] 

25: end while 

26: return n 

27: end method 

28: method Efficient-Find (n) > This speeds up repeated calls. 

29: p <— n 

30: while parent[n] is not None do 

31: n <— parent[n] 

32: end while 

33: while parent [p] ^ n do 

34: p,parent[p] parent[p],n 

35: end while 

36: return n 

37: end method 

38: end class 
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Lines 15 to 19 in NN-CHAIN-CORE clearly find reciprocal nearest neighbors (a, b) in S. 
One important detail is that the index b is preferred in the argmin search in line 16, if the 
minimum is attained at several indices and b realizes the minimum. This can be respected in 
an implementation with no effort, and it ensures that reciprocal nearest neighbors are indeed 
found. That is, the list chain never contains a cycle of length > 2, and a chain — [...,£>, a] 
with reciprocal nearest neighbors at the end will always be extended by b, never with an 
clement c ^ b which coincidentally has the same distance to a. 

After line 19, the chain ends in (6, a, b). The nodes a and b are then joined, and the internal 
variables are updated as usual. 

We now show that the remaining iterations produce the same output as if the algorithm had 
started with the set S' ■— (S \ {a, b}) U {n}, where n is the new node label and the distance 
array d and the size array are updated accordingly. 

The only data which could potentially be corrupted is that the list chain could not contain 
successive nearest neighbors any more, since the new node n could have become the nearest 
neighbor of a node in the list. 

At the beginning of the next iteration, the last elements (b, a, b) are removed from chain. 
The list chain then clearly does not contain a or b at any place any more, since any occurrence 
of a or & in the list would have led to an earlier pair of reciprocal nearest neighbors, before 
(b, a, b) was appended to the list. Hence, chain contains only nodes which really are in S. Let 
e, / be two successive entries in chain, i.e. / is a nearest neighbor of e. Then we know 

d[e, f] < d[e, a] d[a, b] < d[a, e] 

d[e, f] < d[e, b] d[a, b] < d[b, e] 

Together with the reducibility property (1) (for I — a, J = b, K = e), this implies d[e, f] < 
d[e,n]. Hence, / is still the nearest neighbor of e, which proves our assertion. 

We can therefore be sure that the remaining iterations of NN-CHAIN-CORE produce the 
same output as if the algorithm would be run freshly on 5". By the inductive assumption, 
this produces a valid stepwise dendrogram for the set S' with N — 1 nodes. Proposition 3 
carries out the remainder of the proof, as it shows that the first line (a, b, d[a, b]) of the unsortcd 
dendrogram, when it is sorted into the right place in the dendrogram for the nodes in S', is 
a valid stepwise dendrogram for the original set S with N nodes. □ 

Proposition 3. Let (S,d) be a set with dissimilarities (\S\ > I). Fix a distance update 
formula which fulfills the requirements in Theorem 1. Let a,b be two distinct nodes in S 
which are reciprocal nearest neighbors. 

Define S' as (S\ {a, b}) U {n}, where the label n represents the union aUb. Let d' be the up- 
dated dissimilarity matrix for S' , according to the chosen formula. Let L' = ((a^, 6 i; tfj)j=o,...,m) 
be a stepwise dendrogram for S' . Let j be the index such that all Si < d[a,b] for all i < j 
and Si > d[a, b] for all i > j. That is, j is the index where the new row (a, b, d[a, b]) should 
be inserted to preserve the sorting order, giving d[a, b] priority over potentially equal sorting 
keys. Then the array L, which we define as 
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is a stepwise dendrogram for (S,d). 

Proof. Since a and b are reciprocal nearest neighbors at the beginning, the reducibility prop- 
erty (1) guarantees that they stay nearest neighbors after any number of merging steps between 
other reciprocal nearest neighbors. Hence, the first j steps in a dendrogram for S cannot con- 
tain a or b, since these steps all happen at merging dissimilarities smaller than d[a, b]. This is 
the point where we must require that the sorting in line 3 of NN-CHAIN-linkage is stable. 

Moreover, the first j rows of L cannot contain a reference to n: Again by the reducibility 
property, dissimilarities between n and any other node are at least as big as d[a, b]. Therefore, 
the first j rows of L are correct for a dendrogram for S. 

After j steps, we know that no inter-cluster distances in S \ {a, 6} are smaller than d[a, b]. 
Also, d[a, b] is minimal among all distances from a and b, so the row (a,b,d[a,b}) is a valid 
next row in L. 

After this step, we claim that the situation is the same in both settings: The sets S' after 
j steps and the set S after j + 1 steps, including the last one merging a and b into a new 
cluster n, are clearly equal as partitions of the original set. It is required to check that also 
the dissimilarities are the same in both settings. This is where we need the second condition 
in Theorem 1: 

The row (a,b,d[a,b]) on top of the array L' differs from the dendrogram L by j transposi- 
tions, where (a, b, d[a, b]) is moved one step downwards. Each transposition happens between 
two pairs (a, b) and (cii, bi), where all four nodes are distinct, as shown above. The dissimilar- 
ity from a distinct fifth node x to the join a U b does not depend on the merging of a, and bi 
since there is no way in which dissimilarities to aj and bi enter the distance update formula 
Formula (eZ[a, x], d[b, x], d[a, b], size[a], size[b], size[x]). The symmetric statement holds for the 
dissimilarity d[x, a* U bi] . The nodes a,b 7 ai,bi are deleted after the two steps, so dissimilarities 
like d[a, Oj U bi] can be neglected. The only dissimilarity between active nodes which could be 
altered by the transposition is d[aUb, cii U bi]. It is exactly the second condition in Theorem 1 
that this dissimilarity is independent of the order of merging steps. This finishes the proof of 
Theorem 1. □ 

We still have to prove that the requirements of Theorem 1 arc fulfilled by the "single", 
"complete", "average", "weighted" and "Ward" schemes: 

Proof of Proposition 2. It is easy and straightforward to check from the table in Figure 2 
that the distance update schemes in question fulfill the reducibility property. Moreover, the 
table also conveys that the dissimilarities between clusters in the "single", "complete" and 
"average" schemes do not depend on the order of the merging steps. 

For Ward's scheme, the global dissimilarity expression in the third column in Figure 2 
applies only if the dissimilarity matrix consists of Euclidean distances between vectors (which 
is the prevalent setting for Ward's method). For a general argument, note that the global 
cluster dissimilarity for Ward's method can also be expressed by a slightly more complicated 
expression: 



d(A, B) = 1- U £ £ d(a, by M £ £ d(a, a') 2 ~^EE d ^ "A 

\ 1 1 1 1 \ aeAb&B 1 1 aeAa'EA 1 1 btBb'tB / 

This formula can be proved inductively from the recursive distance update formula for Ward's 
method, hence it holds independently of whether the data is Euclidean or not. This proves 
that the dissimilarities in Ward's scheme arc also independent of the order of merging steps. 

Dissimilarities in the "weighted" scheme, however, do in general depend on the order of 
merging steps. However, the dissimilarity between joined nodes I Li J and K U L is always the 
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Figure 6 The single linkage algorithm. 



procedure MST-linkage(S, d) > S: node labels, d: pairwisc dissimilarities 

L <- MST-linkage-core(S, d) 
Stably sort L with respect to the third column. 

L <— Label(L) > Find node labels from cluster representatives, 

return L 
end procedure 

procedure MST- LINKAGE-CO Re(jSo, d) > S a : node labels, d: pairwise dissimilarities 
L<-[] 

c «— (any element of So) > c: current node 

D [s] «- oo for seSo\ {c} 
for i in (1, ... , |5o| — 1) do 

5, <- \ {c} 

for s in Sj do 

Di[s] <- min{Di-i[s],d[s,c]} 

end for 

n <— argmin sgS . Z?j[s] > n: new node 

Append (c, n, D^n]) to L 
c n 
end for 

return L > an unsorted dendrogram 

end procedure 



mean dissimilarity if] + d[I, L] + d{J, K] + d[J, L]), independent of the order of steps, 

and this is all that is required for Proposition 2. □ 

3.3 The single linkage algorithm 

In this section, we present and prove correctness of a fast algorithm for single linkage clus- 
tering. Gowcr and Ross (1969) observed that a single linkage dendrogram can be obtained 
from a minimum spanning tree (MST) of the weighted graph which is given by the complete 
graph on the singleton set S with the dissimilarities as edge weights. The algorithm here 
was originally described by Rohlf (1973) and is based on Prim's algorithm for the MST (see 
Cormcn et ah, 2009, §23.2). 

The single linkage algorithm MST-linkage is given in Figure 6. In the same way as the NN- 
chain algorithm, it consists of a core algorithm MST-linkage-CORE and two postprocessing 
steps. The output structure of the core algorithm is again an unsorted list of clustering 
steps with node representatives instead of unique labels. As will be proved, exactly the same 
postprocessing steps can be used as for the NN-chain algorithm. 

Rohlf's algorithm in its original version is a full Prim's algorithm and maintains enough 
data to generate the MST. He also mentions a possible simplification which does not do 
enough bookkeeping to generate an MST but enough for single linkage clustering. It is this 
simplification that is discussed in this paper. We prove the correctness of this algorithm for 
two reasons: 

• Since the algorithm MST-linkage-CORE does not generate enough information to re- 
construct a minimum spanning tree, one cannot refer to the short proof of Prim's algo- 
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rithm in any easy way to establish the correctness of MST-linkage. 

• Like for the NN-chain algorithm in the last section, it is not clear a priori that the 
algorithm resolves ties correctly. A third algorithm can serve as a warning here (see 
Section 5 for more details): There is an other fast algorithm for single linkage cluster- 
ing, Sibson's SLINK algorithm (Sibson, 1973). More or less by coincidence, all three 
algorithms NN-CHAIN-CORE, MST-linkage-CORE and SLINK generate output which 
can be processed by exactly the same two steps: sorting followed by Label. In case 
of the SLINK algorithm this works fine if all dissimilarities are distinct but produces 
wrong stepwise dendrograms in situations when two merging dissimilarities are equal. 
There is nothing wrong with the SLINK algorithm, however. Sibson supplied a proof for 
the SLINK algorithm in his paper (Sibson, 1973), but it is written for a (non-stepwise) 
dendrogram as the output structure, not for a stepwise dendrogram. Hence, the addi- 
tional information which is contained in a stepwise dendrogram in the case of ties is not 
provided by all, otherwise correct algorithms. 

This should be taken as a warning that ties demand more from an algorithm and must be 
explicitly taken into account when we prove the correctness of the MST-linkage algorithm 
below. 

Theorem 4. The algorithm MST-LINKAGE yields an output which can also be generated by 
Primitive_clustering . 

We do not explicitly refer to Prim's algorithm in the following, and we make the proof self- 
contained, since the algorithm does not collect enough information to construct a minimum 
spanning tree. There are unmistakable similarities, of course, and the author got most of the 
ideas for this proof from Prim's algorithm (see Cormen ct al., 2009, §23.2). 

Let us first make two observations about the algorithm MST-linkage. 

(a) Starting with the full initial set So, the algorithm MST-linkage-CORE chooses a "cur- 
rent node" c in each step and removes it from the current set Si in every iteration. Let 
Sf := S \ Si be the complement of the current set Si. Then Di[s] (s e Si) is the 
distance from Sf to s, i.e. 

Di[s] = mm d[s, t]. 

(b) Let L be the output of MST-LINKAGE-CORE^, d). The 2i entries in the first two 
columns and the first i rows contain only i + 1 distinct elements of S, since the second 
entry in one row is the first entry in the next row. 

We prove Theorem 4 by the following stronger variant: 

Theorem 5. Let L be the output of MST-linkage-CORe(S i , d). For all n < \S\, the first n 
rows of L are an unsorted single linkage dendrogram for the n+ 1 points of S in this list (see 
Observation (b)). 

Proof. We proceed by induction. After the first iteration, the list L contains one triple 
(do, bo, So). So = Di[bo] is clearly the dissimilarity d[ao,&o], since the array D\ contains the 
dissimilarities to ao after the first iteration (Observation (a)). 

Let (ao, bo, So), ■■■ , (a n , b n , S n ) be the first n + 1 rows of L. We sort the rows with a stable 
sorting algorithm as specified in MST-linkage. We leave the postprocessing step Label out 
of our scope and work with the representatives a i; 6j for the rest of the proof. 
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Let s(0), . . . , s(ri) be the stably sorted indices (i.e. 8 S ^) < <5 s (i+i) for all i and s(z) < s(z+ 1) 
if ^s(i) = <^s(i+i))- Let k be the sorted index of the last row n. Altogether, we have a sorted 
matrix 





a s(0) 


b s (o) 


<>s(0) 




a s (k-i) 


b s (k-i) 


^s(fe-l) 




a s (k) 


b s (k) 


^s(fe) 




a s(k + l) 


b s (k+i) 


8 s (k+i) 






b s ( n ) 





The new row is at the index k, i.e. (a s (fc), 6 s (fe), # s (fc)) = (a n ,b n , 8 n ). The matrix without 
the fc-th row is a valid stepwise, single linkage dendrogram for the points a®, . . . ,a n , by the 
induction hypothesis. (Recall that bi = Oj+i.) Our goal is to show that the matrix with its 
fc-th row inserted yields a valid single linkage dendrogram on the points ao, . . . ,a n ,b n . 

First step: rows to k — 1. The distance 8 n is the minimal distance from b n to any 
of the points ao,.. . ,a n . Therefore, the dendrograms for the sets S~ := {a , . . . ,a n } and 
S + := S~ U {b n } have the same first k steps, when all the inter-cluster distances are smaller 
than or equal to S n . (If the distance 8 n occurs more than once, i.e. when 8 s ^-i) = we 
assume by stable sorting that the node pairs which do not contain b n are chosen first.) 

Therefore, the first k rows are a possible output of Primitive_CLUSTERING in the first k 
steps. After this step, we have the same partial clusters in S + as in the smaller data set, plus 
a singleton {&„}. 

Second step: row k. The distance 8 n from b n to some point ao,...,a n is clearly the 
smallest inter-cluster distance at this point, since all other inter-cluster distances are at least 
^s(fc+i)i which is greater than 8 S ^) = 8 n . Since the output row is (a n ,b n ,8 n ), it remains to 
check that the distance S n is realized as the distance from b n to a point in the cluster of a n , 
i.e. that a n is in a cluster with distance 8 n to b n . 

The clusters mentioned in the last sentence refer to the partition of S + which is generated 



by the relations a s ( ) ~ 6 S ( ), . . . ,a s (fe-i) ~ ^s(fe-i)- Since we have bi - 
S + consists of contiguous chains in the original order of points ao, ai, 
The diagram below visualizes a possible partition after k steps. 



aj+i, the partition of 

• • ; ^n; b n . 



In this particular example, the distances So, 5i and Sn are among the first k smallest, while 
82, £3 and S 5 come later in the sorted order. 

Let 8 n be realized as the distance between b n and a m for some m < n. Then the dis- 
similarities between consecutive points in the sequence a m , b m = a m+ \, b m+1 = a m+2 , . . . , 
6„_i = a n must be less than or equal to 8 n ; otherwise b n and the dissimilarity 8 n = d[b n ,a m ] 
would have been chosen first over these other dissimilarities in one of the first k steps. Since 
the dissimilarities of all pairs (a,i,bi) in this chain are not more than 8 n , they are contained 
in the first k sorted triples. Hence, a m and a n have been joined into a cluster in the first k 
steps, and a n is a valid representative of the cluster that also contains a m . 

Note that the argument in the last paragraph is the point where we need that the sorting 
in line 3 of MST-linkage is stable. Otherwise it could not be guaranteed that a m and a n 
have been joined into a cluster before b n is added. 

Third step: rows k+ 1 to n. Here is the situation after row k: We have the same clusters 
in S + as after k steps in the smaller data set S, except that the last cluster (the one which 
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contains a n ) additionally contains the point b n . In a diagram: 

( qq «i a2] [a3]( a4 05] ••• [■■■ a m ■■■ a n K 

The inter-cluster distances in S + from the cluster with b n to the other clusters might be 
smaller than without the point b n in S~ . We show, however, that this does not affect the 
remaining clustering steps: 

In each step r > k, we have the following situation for some x < y < s(k). The point b n 
might or might not be in the same cluster as b s r r y 



l s(r) 



t>s(r) — a s{r) + l b n -l 



Let the distance <$ s ( r ) be realized as the distance from b s t r -\ to a y . From Observation (a) and 
line 10 in MST-linkage-CORE, we know that this distance is minimal among the distances 
from X := {a , . . . , a s ( r )} to all other points in S \ X. In particular, the distance from X to 
b n G So \ X is not smaller than 6 s ( r y 

This proves that the addition of b n in step k does not change the single linkage clustering 
in any later step r > k. This completes the inductive proof of Theorem 5 □ 



4 Performance 

In this section, we compare the performance of the algorithms and give recommendations on 
which algorithm to choose for which clustering method. We compare both the theoretical, 
asymptotic worst-case performance, and the use-case performance on a range of synthetic 
random data sets. 



4.1 Asymptotic worst-case performance 

Let N denote the problem size, which is in this case the number of input data points. The 
input size is (") e Q(N 2 ). 

The asymptotic run-time complexity of MST-linkage-CORE is obviously 6(iV 2 ), since 
there are two nested levels of loops in the algorithm (line 8 and implicitly line 10). The 
run-time complexity of the NN-CHAIN-CORE algorithm is also <d(N 2 ) (Murtagh, 1985, page 
86). Postprocessing is the same for both algorithms and is less complex, namely 0(N log N) 
for sorting and Q(N) for Label, so the overall complexity is 0(iV 2 ). This is optimal (in the 
asymptotic sense): the lower bound is also quadratic since all Q(N 2 ) input values must be 
processed. 

The NN-chain algorithm needs a writable working copy of the input array to store inter- 
mediate dissimilarities and otherwise only Q(N) additional memory. 

The generic algorithm has a best-case time complexity of 6(iV 2 ), but without deeper anal- 
ysis, the worst-case complexity is 0(N 3 ). The bottleneck is line 15 in Generic linkage: 

In O(N) iterations, this line might be executed up to O(N) times and does a minimum search 
over O(N) elements, which gives a total upper bound of 0(N 3 ). This applies for all clustering 
schemes except single linkage, where the loop starting at line 14 is never executed and thus 
the worst-case performance is 0(-/V 2 ). The memory requirements for the generic algorithm are 
similar to the NN-chain algorithm: a working copy of the dissimilarity array and additionally 
only 0(iV) temporary memory. 

In contrast, the MST algorithm does not write to the input array d. All other temporary 
variables are of size O(N). Hence, MST-linkage requires no working copy of the input 
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array and hence only half as much memory as Generic_LINKAGE and NN-CHAIN-linkage 
asymptotically. 

Andcrbcrg's algorithm (Anderberg, 1973, pages 135-136) has the same asymptotic bounds 
as our generic algorithm. The performance bottleneck are again the repeated minimum 
searches among the updated dissimilarities. Since the generic algorithm defers minimum 
searches to a later point in the algorithm (if they need to be performed at all, by then), 
there are at least as many minimum searches among at least as many elements in Ander- 
berg's algorithm as in the generic algorithm. The only point where the generic algorithm 
could be slower is the maintenance of the priority queue with nearest neighbor candidates, 
since this does not exist in Anderberg's algorithm. The bottleneck here are potentially up to 

0(N 2 ) updates of a queue in line 36 of Generic LINKAGE. In the implementation in the 

next section, the queue is realized by a binary heap, so an update takes 0(log N) time. This 
could potentially amount to O (TV 2 log N) operations for maintenance of the priority queue. 
However, a reasonable estimate is that the saved minimum searches in most cases save more 
time than the maintenance of the queue with O(N) elements costs, and hence there is a good 
reason to believe that the generic algorithm is at least as fast as Anderberg's algorithm. 

Note that the maintenance effort of the priority queue can be easily reduced to 0(N 2 ) 
instead of 0(N 2 log N) worst case: 

• A different priority queue structure can be chosen, where the "decrease-key" operation 
takes only O(l) time. (Note that the bottleneck operation in line 36 of Generic LINK- 
AGE never increases the nearest-neighbor distance, only decreases it.) The author did 
not test a different structure since a binary heap convinces by its simple implementation. 

• Changed keys (minimal distances) need not be updated in the priority queue imme- 
diately. Instead, the entire queue might be resorted/regenerated at the beginning of 
every iteration. This takes N — 1 times O(N) time with a binary heap. Although this 
lowers the theoretical complexity for the maintenance of the binary queue, it effectively 
slowed down the algorithms in practice by a small margin. The reason is, of course, that 
the number and complexity of updates of the priority queue did by far not reach their 
theoretical upper bound in our test data sets (see below). Altogether, the maintenance 
of the priority queue, as proposed in Figure 3 seems quite optimal from the practical 
perspective. 

4.2 Use-case performance 

In addition to the theoretical, asymptotic and worst-case considerations, we also measured 
the practical performance of the algorithms. Figure 7 shows the run-time of the algorithms 
for a number of synthetic test data sets (for details see below) . The solid lines are the average 
over the data sets. (The graphs labeled "Day-Edelsbrunner" are discussed in Section 5.) The 
lightly colored bands show the range from minimum to maximum time over all data sets for 
a given number of points. 

The following observations can be made: 

• For single linkage clustering, the MST-algorithm is clearly the fastest one. Together 
with the fact that it has only half the memory requirements of the other algorithms (if 
the input array is to be preserved), and thus allows the processing of larger data sets, 
the MST-algorithm is clearly the best choice for single linkage clustering. 

• For the clustering schemes without inversions (all except "centroid" and "median"), the 
generic algorithm, the NN-chain algorithm and Anderberg's algorithm have very similar 
performance. 
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Method: "single" 



Method: "average" 

"complete", "weighted" and 
"Ward" look very similar. 



Method: "centroid" 
"median" looks very similar. 




Number of points 

Figure 7: Performance of several SAHN clustering algorithms. Legend: ■ Generic algorithm 
(Figure 3), ■ Anderberg (1973, pages 135-136), ■ NN-chain algorithm (Figure 4), 
■ MST-algorithm (Figure 6), ■ Day and Edelsbrunner (1984, Table 5). 
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The NN-chain algorithm is the only one with guaranteed 0(N 2 ) performance here. We 
can conclude that the good worst-case performance can be had here without any cut- 
backs to the use-case performance. 

• For the "centroid" and "median" methods, we see a very clear disadvantage to An- 
derberg's algorithm. Here, the worst case cubic time complexity occurs already in the 
random test data sets. This happens with great regularity over the full range of input 

sizes. Our Generic linkage algorithm, on the other hand, does not suffer from this 

weakness: Even though the theoretical worst-case bounds are the same, the complexity 
does not raise above the quadratic behavior in our range of test data sets. Hence, we 
have grounds to assume that Generic_linkage is much faster in practice. 

4.3 Conclusions 

Based on the theoretical considerations and use-case tests, we can therefore recommend algo- 
rithms for the various distance update schemes as follows: 

• "single" linkage clustering: The MST-algorithm is the best, with respect to worst-case 
complexity, use-case performance and memory requirements. 

• "complete", "average", "weighted", "ward": The NN-chain algorithm is preferred, since 
it guarantees 0(N 2 ) worst case complexity without any disadvantage to practical per- 
formance and memory requirements. 

• "centroid", "median": The generic clustering algorithm is the best choice, since it can 
handle inversions in the dendrogram and the performance exhibits quadratic complexity 
in all observed cases. 

Of course, the timings in the use-case tests depend on implementation, compiler optimiza- 
tions, machine architecture and the choice of data sets. Nevertheless, the differences between 
the algorithms are very clear here, and the comparison was performed with careful implemen- 
tations in the identical environment. 

The test setup was as follows: All algorithms were implemented in CH — h with an interface 
to Python (van Rossum et al.) and the scientific computing package NumPy (Num) to handle 
the input and output of arrays. The test data sets are samples from mixtures of multivariate 
Gaussian distributions with unity covariance matrix in various dimensions (2, 3, 10, 200) with 
various numbers of modes (1, 5, [%/iV]), ranging from N = 10 upwards until memory was 
exhausted (N = 20000 except for single linkage). The centers of the Gaussian distributions 
are also distributed by a Gaussian distribution. Moreover, for the methods for which it makes 
sense (single, complete, average, weighted: the "combinatorial" methods), we also generated 
10 test sets per number of input points with a uniform distribution of dissimilarities. 

The timings were obtained on a PC with an Intel dual-core CPU T7500 with 2.2 GHz clock 
speed and 4GB of RAM and no swap space. The operating system was Ubuntu 11.04 64-bit, 
Python version: 2.7.1, NumPy version: 1.5.1, compiler: GNU C++ compiler, version 4.5.2. 
Only one core of the two available CPU cores was used in all computations. 

5 Alternative algorithms 

The MST algorithm has the key features that it (1) needs no working copy of the 0(iV 2 ) input 
array and only O(iV) working memory (2) is fast since it reads every input dissimilarity only 
once and otherwise deals only with 6(iV) memory. There is a second algorithm with these 
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characteristics, Sibson's SLINK algorithm (Sibson, 1973). It is based on the insight that a 
single linkage dendrogram for N + 1 points can be computed from the dendrogram of the first 
N points plus a single row of distances (d[N, 0], . . . , d[N, N — 1]). In this fashion, the SLINK 
algorithm even reads the input dissimilarities in a fixed order, which can be an advantage 
over the MST algorithm if the favorable input order can be realized in an application, or if 
dissimilarities do not fit into random-access memory and are read from disk. 

However, there is one important difference: even though the output data format looks 
deceptively similar to the MST algorithm (the output can be converted to a stepwise den- 
drogram by exactly the same process: sorting with respect to dissimilarities and a union-find 
procedure to generate node labels from cluster representatives), the SLINK algorithm can- 
not handle ties. This is definite, since e.g. the output in the example situation on page 5 is 
the same in all three cases, and hence no postprocessing can recover the different stepwise 
dendrograms. 

There is an easy way out by specifying a secondary order 



to make all dissimilarities artificially distinct. In terms of performance, the extra comparisons 
put a slight disadvantage on the SLINK algorithm, according to the author's experiments. 
However, the difference is not much, and the effect on timings may be compensated or even 
reversed in a different software environment or when the input order of dissimilarities is in 
favor of SLINK. Hence, the SLINK algorithm is a perfectly fine tool, as long as care is taken 
to make all dissimilarities unique. 

The same idea of generating a dendrogram inductively is the basis of an algorithm by 
Dcfays (1977). This paper is mostly cited as a fast algorithm for complete linkage clustering. 
However, it definitely is not an algorithm for complete linkage clustering, as the complete 
linkage method is commonly defined, in this paper and identically elsewhere. 

An algorithm which is interesting from the theoretical point of view is given by Day and 
Edelsbrunner (1984, Table 5). It uses N priority queues for the nearest neighbor of each 
point. By doing so, the authors achieve a worst-case time complexity of 0(N 2 log N), which 
is better than the existing bound 0(N 3 ) for the schemes where the NN-chain algorithm 
cannot be applied. The overhead for maintaining a priority queue for each point, however, 
slows the algorithm down in practice. The performance measurements in Figure 7 include 
the Day-Edelsbrunner algorithm. Day and Edelsbrunner write their algorithm in general 
terms, for any choice of priority queue structure. We implemented the algorithm for the 
measurements in this paper with binary heaps, since these have a fixed structure and thus 
require the least additional memory. But even so, the priority queues need additional memory 
of order Q(N 2 ) for their bookkeeping, which can also be seen in the graphs since they stop 
at fewer points, within the given memory size of the test. The graphs show that even if the 
Day-Edelsbrunner algorithm gives the currently best asymptotic worst-case bound for the 
"centroid" and "median" methods, it is inefficient for practical purposes. 

Kfivanek (1990, §11) suggested to put all (^) dissimilarity values into an (a, 6)-tree data 
structure. He claims that this enables hierarchical clustering to be implemented in 0(N 2 ) 
time. Kfivanek's conceptually very simple algorithm relies on the fact that m insertions into 
an (a, &)-tree can be done in O(m) amortized time. This is only true when the positions, 
where the elements should be inserted into the tree, are known. Searching for these positions 
takes O(logiV) time per element, however. (See Mehlhorn and Tsakalidis (1990, §2.1.2) for 
an accessible discussion of amortized complexity for (2,4)-trees; Huddlcston and Mehlhorn 
(1982) introduce and discuss (a, &)-trees in general.) Kfivanek did not give any details of his 
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analysis, but based on his short remarks, the author cannot see how Kfivanek's algorithm 
achieves 0(N 2 ) worst-case performance for SAHN clustering. 

6 Extension to vector data 

If the input to a SAHN clustering algorithm is not the array of pairwise dissimilarities but 
N points in a D-dimcnsional real vector space, the lower bound il(N 2 ) on time complexity 
does not hold any more. Since much of the time in an SAHN clustering scheme is spent on 
nearest-neighbor searches, algorithms and data structures for fast nearest-neighbor searches 
can potentially be useful. The situation is not trivial, however, since (1) in the "combinatorial" 
methods (e.g. single, complete, average, weighted linkage) the inter-cluster distances are not 
simply defined as distances between special points like cluster centers, and (2) even in the 
"geometric" methods (the Ward, centroid and median schemes), points are removed and 
new centers added with the same frequency as pairs of closest points are searched, so a 
dynamic nearest-neighbor algorithm is needed, which handles the removal and insertion of 
points efficiently 

Moreover, all known fast nearest-neighbor algorithms lose their advantage over exhaustive 
search with increasing dimensionality. Additionally, algorithms will likely work for one metric 
in WL D but not universally Since this paper is concerned with the general situation, we do 
not go further into the analysis of the "stored data approach" (Anderberg, 1973, §6.3). We 
only list at this point what can be achieved with the algorithms from this paper. This will 
likely be the best solution for high-dimensional data or general-purpose algorithms, but there 
are better solutions for low-dimensional data outside the scope of this paper. The suggestions 
below are at least helpful to process large data sets since memory requirements are of class 
Q(ND), but they do not overcome their Q(N 2 ) lower bound on time complexity. 

• The MST algorithm for single linkage can compute distances on-the-fly. Since every 
pairwise dissimilarity is read in only once, there is no performance penalty compared to 
first computing the whole dissimilarity matrix and then applying the MST algorithm. 
Quite the contrary, computing pairwise distances in-process can result in faster execu- 
tion since much less memory must be reserved and accessed. The MST algorithm is 
suitable for any dissimilarity measure which can be computed from vector representa- 
tions (that is, all scale types are possible, e.g. M-valued measurements, binary sequences 
and categorical data). 

• The NN-chain algorithm is suitable for the "Ward" scheme, since inter-cluster distances 
can be defined by means of centroids as in Figure 2. The initial inter-point dissimilarities 
must be Euclidean distances (which is anyway the only setting in which Ward linkage 
describes a meaningful procedure). 

• The generic algorithm is suitable for the "Ward", "centroid" and "median" scheme on 
Euclidean data. There is a simpler variant of the Generic linkage algorithm in Sec- 
tion 6, which works even faster in this setting. The principle of the algorithm Generic 

linkage variant is the same: each array entry mindist[x] maintains a lower bound on 

all dissimilarities d[x, y] for nodes with label y > x. The Generic linkage algorithm 

is designed to work efficiently with a large array of pairwise dissimilarities. For this pur- 
pose, the join of two nodes a and b re-uses the label b, which facilitates in-place updating 
of the dissimilarity array in an implementation. The Generic linkage_VARIANT al- 
gorithm, in contrast, generates a unique new label for each new node, which is smaller 
than all existing labels. Since the new label is at the beginning of the (ordered) list of 
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Figure 8 The generic clustering algorithm (variant). 

l: procedure Generic linkage variant(./V, d) 

d is either an array or a function which computes dissimilarities from cluster centers. 

(Lines 2 to 13 arc the same as in Generic_linkage.) 
5: for x in S \ {N - 1} do 

14: while b S do 

: > Recalculation of nearest neighbors, if necessary. 

21: end while 

22: Remove a and b from Q. 

23: Append (a, b, S) to L. 

24: Create a new label n < i 

25: size[n] size [a] + size[b] 

26: S 4- (S \ {a, b}) U {n} 

27: for x in S \ {n} do > Extend the distance information. 

28: d[x, n] 4— rf[n,a;] <— FoRMULA(c£[a, a;], d[b, x], d[a, b], size[a], size[b], size[x\) 

29: end for 
or 

27: Compute the cluster center for n as in Figure 2. 

30: n_nghbr [n] argmin 2 . >n d[n, x] 

31: Insert (n, d[n, n_nghbr[n]]) into mindist and Q 

32: end for 

33: return L 

34: end procedure 



nodes and not somewhere in the middle, the bookkeeping of nearest neighbor candidates 

and minimal distances is simpler in Generic linkage_variant: in particular, the 

two loops in lines 28-38 of Generic linkage can be disposed of entirely. Moreover, 

experiments show that Generic linkage variant needs much less recalculations of 

nearest neighbors in some data sets. However, both algorithms are similar, and which 
one is faster in an implementation seems to depend strongly on the actual data structures 
and their memory layout. 

Another issue which is not in the focus of this paper is that of parallel algorithms. For the 
"stored matrix approach", this has a good reason since the balance of memory requirements 
versus computational complexity does not make it seem worthwhile to attempt parallelization 
with current hardware. This changes for vector data, when the available memory is not the 
limiting factor and the run-time is pushed up by bigger data sets. In high-dimensional vector 
spaces, the advanced clustering algorithms in this paper require little time compared to the 
computation of inter-cluster distances. Hence, parallelizing the nearest- neighbor searches with 
their inherent distance computations appears a fruitful and easy way of sharing the workload. 
The situation becomes less clear for low-dimensional data, however. 
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7 Conclusion 



Among the algorithms for sequential, agglomcrativc, hierarchic, nonovcrlapping (SAHN) 
clustering on data with a dissimilarity index, three current algorithms are most efficient: 
Rohlf's algorithm MST-linkage for single linkage clustering, Murtagh's algorithm NN- 
CHAIN-linkage for the "complete", "average", "weighted" and "Ward" schemes, and the 

author's Generic linkage algorithm for the "centroid" and "median" schemes and the 

"flexible" family. The last algorithm can also be used for an arbitrary distance update for- 
mula. There is even a simpler variant Generic LINKAGE variant, which seems to require 

less internal calculations, while the original algorithm is optimized for in-placc updating of 

a dissimilarity array as input. The Generic LINKAGE algorithm and its variant are new; 

the other two algorithms were described before, but for the first time they are proved to be 
correct. 
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